
setwd(pathDATA)

coeff <- summary(sp.weib.v$mle)[[6]][,1]
betas <- coeff[1:11]
rho <- coeff[length(summary(sp.weib.v$mle)[[6]][,1])-1]
shape <- coeff[length(summary(sp.weib.v$mle)[[6]][,1])]
data <- apply(sp.weib.v$data[,3:12],2,FUN=median)
data <- rbind(data,data,data,data)
W <- matrix(1/3,nrow=4,ncol=4)
diag(W) <- 0

eye<-matrix(nrow=4, ncol=4, 0)
diag(eye) <- 1


weibullsim<-function(betas,data,rho,W,shape,eye){
	constant	<- 1
	betas		<- as.numeric(betas)
	X			<- as.matrix(data.frame(constant,data))
	gumbel		<- log(rweibull(dim(X)[1], shape=1, scale = 1))
    error		<- (shape)*gumbel
	duration	<- exp(solve(eye-(rho*W))%*%(X%*%betas+error))
	return(duration)
	}
	


results.rho <-matrix(0,ncol=10000,nrow=4)
results.rho0 <-matrix(0,ncol=10000,nrow=4)
results.rhominus1 <-matrix(0,ncol=10000,nrow=4)
results.rho1 <-matrix(0,ncol=10000,nrow=4)
	
for(i in 1:10000){
	results.rho[,i] <- weibullsim(betas=betas,data=data,rho=rho,W=W,shape=shape,eye=eye)
	results.rho0[,i] <- weibullsim(betas=betas,data=data,rho=0,W=W,shape=shape,eye=eye)
	results.rhominus1[,i] <- weibullsim(betas=betas,data=data,rho=-0.1,W=W,shape=shape,eye=eye)
	results.rho1[,i] <- weibullsim(betas=betas,data=data,rho=0.1,W=W,shape=shape,eye=eye)
	}


#Figure 5
pdf(file="rho0.pdf",width=5,height=5)
hist(results.rho0,freq=FALSE,ylim=c(0,0.01),xlim=c(0,1000),breaks=100,col="grey",border="grey",ann=FALSE)
lines(c(median(results.rho0),median(results.rho0)),c(0,-1))
title(ylab="Density",xlab="Duration in days")
text(c(median(results.rho0)),0,lab="Median",pos=3)
dev.off()

pdf(file="rho1.pdf",width=5,height=5)
hist(results.rho1,freq=FALSE,ylim=c(0,0.01),xlim=c(0,1000),breaks=100,col="grey",border="grey",ann=FALSE)
lines(c(median(results.rho1),median(results.rho1)),c(0,-1))
title(ylab="Density",xlab="Duration in days")
text(c(median(results.rho1)),0,lab="Median",pos=3)
dev.off()


pdf(file="rhominus1.pdf",width=5,height=5)
hist(results.rhominus1,freq=FALSE,ylim=c(0,0.01),xlim=c(0,1000),breaks=100,col="grey",border="grey",ann=FALSE)
lines(c(median(results.rhominus1),median(results.rhominus1)),c(0,-1))
title(ylab="Density",xlab="Duration in days")
text(c(median(results.rhominus1)),0,lab="Median",pos=3)
dev.off()

median(results.rhominus1)
median(results.rho0)
median(results.rho1)


#hist(results.rho2,freq=FALSE,ylim=c(0,0.002),xlim=c(0,10000),breaks=100)

	
